% Copyright (C) 2009,2010,2011,2012  Marco Restelli
%
% This file is part of:
%   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
%
% LDGH is free software: you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% LDGH is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
% License for more details.
%
% You should have received a copy of the GNU General Public License
% along with LDGH. If not, see <http://www.gnu.org/licenses/>.
%
% author: Marco Restelli                   <marco.restelli@gmail.com>


function grid_oct2gmsh(p,e,t,gridfile)
% FIX THE COMMENTS
% [p,e,t] = grid_gmsh2oct(d,gridfile,renumber)
%
% This function reads a .msh gmsh grid from file "gridfile" and
% reformats it in the usual [p,e,t] "Matlab-like" format.
%
% For 2D meshes, the third coordinate is ignored.
%
% renumber is optional and, if present, it defines a renumbering for
% the boundary labels in e. The length of renumber must be (at least)
% equal to the maximum boundary label used in gridfile. Consider that:
% a) setting renumber(i)<0 the sides with boundary label i are
%   discarded (this is useful, for instance, to drop internal
%   boundaries)
% b) if it's known that an intermediate label i is not used, the
%   corresponding element renumber(i) should be NaN, in order to
%   detect any error.
%
% LIMITATIONS: 
% - only one label for the domains is recognized
% - the curvilinear abscissa in e is not set

 fid = fopen(gridfile,'w');

 % Header
 line = '$MeshFormat';
 fprintf(fid,'%s\n',line);
 line = '2.2 0 8';
 fprintf(fid,'%s\n',line);
 line = '$EndMeshFormat';
 fprintf(fid,'%s\n',line);

 % Physical Names
 line = '$PhysicalNames';
 fprintf(fid,'%s\n',line);
 fprintf(fid,'%1d\n',2);
 fprintf(fid,'%1d %1d %s\n',2,1,'"Dirichlet"');
 fprintf(fid,'%1d %1d %s\n',3,2,'"Regione_intera"');
 line = '$EndPhysicalNames';
 fprintf(fid,'%s\n',line);

 % Vertices
 line = '$Nodes';
 fprintf(fid,'%s\n',line);
 fprintf(fid,'%1d\n',size(p,2));
 for i=1:size(p,2)
   fprintf(fid,'%1d %.15e %.15e %.15e\n',i,p(1,i),p(2,i),p(3,i));
 end
 line = '$EndNodes';
 fprintf(fid,'%s\n',line);

 % Boundary faces and elements
 line = '$Elements';
 fprintf(fid,'%s\n',line);
 fprintf(fid,'%1d\n',size(e,2)+size(t,2));
 for i=1:size(e,2)
   fprintf(fid,'%1d %1d %1d %1d %1d %1d %1d\n',i,2,1,1, ...
           e(1,i),e(2,i),e(3,i));
 end
 for i=1:size(t,2)
   ii = size(e,2) + i;
   fprintf(fid,'%1d %1d %1d %1d %1d %1d %1d %1d\n',ii,4,1,2, ...
           t(1,i),t(2,i),t(3,i),t(4,i));
 end
 line = '$EndElements';
 fprintf(fid,'%s\n',line);

 fclose(fid);

return

